{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further Hypothesis Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select this cell and type Ctrl-Enter to execute the code below.\n",
    "\n",
    "library(tidyverse)\n",
    "\n",
    "set_plot_dimensions <- function(width_choice, height_choice) {\n",
    "    options(repr.plot.width=width_choice, repr.plot.height=height_choice)\n",
    "}\n",
    "\n",
    "cbPal <- c(\"#E69F00\", \"#56B4E9\", \"#009E73\", \"#F0E442\", \"#CC79A7\", \"#0072B2\", \"#D55E00\")\n",
    "\n",
    "set_plot_dimensions(5, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You should see \"Attaching packages\" and some ticks by the packages loaded.\n",
    "# The \"Conflicts\" aren't a problem.\n",
    "\n",
    "# Other problems loading the library? Try running this cell.\n",
    "\n",
    "install.packages('tidyverse')\n",
    "\n",
    "library(tidyverse)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8 - Goodness of fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run this cell to load the data.\n",
    "\n",
    "data <- read_csv(\"stars.csv\")\n",
    "\n",
    "type_key <- c('Brown Dwarf', 'Red Dwarf', 'White Dwarf', 'Main Sequence', 'Supergiant','Hypergiant')\n",
    "spectral_classes <- c('O','B','A','F','G','K','M')\n",
    "\n",
    "data$type <- factor(data$type)\n",
    "data$spectral_class <- factor(data$spectral_class, levels=spectral_classes)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your colleague Althea is not a fan of Prof. Xu's temperature scheme. She has her own ideas about how the star classification should be revised.\n",
    "\n",
    "According to Althea's new theory of stellar evolution, there are two subtypes of hypergiant star, with very different temperature distributions. \n",
    "\n",
    "She says that our galaxy has approximately equal numbers of each subtype. \n",
    "\n",
    "If this theory is true, hypergiant temperatures should be distributed according to:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "P(3000K \\leq \\text{T} \\lt 4000K) = 0.5 \\\\\n",
    "P(4000K \\leq \\text{T} \\lt 40000K) = 0.5\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with uniform temperature distributions within each of the two subtypes.\n",
    "\n",
    "Althea asks you to check whether her theory agrees with your data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: do hypergiants in the observed data set fit Althea's theory?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pdf associated with the theory is piecewise uniform. We will start by constructing it as a bespoke function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "althea_pdf_single <- function(x){\n",
    "    if(x<3000) return(0)\n",
    "    else if(x<4000) return(0.5/(4000-3000))\n",
    "    else if(x<40000) return(0.5/(40000-4000))\n",
    "    else return(0)\n",
    "}\n",
    "\n",
    "# We can use the Vectorize() function to make a vector-friendly version of this pdf:\n",
    "althea_pdf <- Vectorize(althea_pdf_single)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "x <- seq(0,50000)\n",
    "y <- althea_pdf(x)\n",
    "\n",
    "xlab <- 'temperature'\n",
    "ylab <- 'probability density'    \n",
    "\n",
    "plot(x, y, type='l', col='orange', xlab=xlab, ylab=ylab, main='')\n",
    "legend('topright', legend='theoretical pdf', \n",
    "       col='orange', lty=1, cex=0.8, y.intersp=2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first glance, the histogram for the hypergiants (type 5) does appear to have a similar shape to this pdf:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "observed <- filter(data, type == 5)$temperature\n",
    "n <- length(observed)\n",
    "\n",
    "xlab <- 'temperature'\n",
    "ylab <- 'freq'\n",
    "\n",
    "bins <- seq(0, 50000, 500)\n",
    "\n",
    "hist(observed, breaks=bins, xlab=xlab, ylab=ylab, col=cbPal[6], main='')\n",
    "legend('topright', legend='observed data', fill=cbPal[6], cex=0.8, y.intersp=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see more clealy, we can overlay a random sample of the same size ($n$=40), drawn from the theoretical distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r <- runif(n=n,min=0,max=1)\n",
    "sample <- rep(NA,n)\n",
    "sample[r<0.5] <- 3000 + (4000-3000)*(r[r<0.5]/0.5)\n",
    "sample[r>=0.5] <- 4000 + (40000-4000)*((r[r>=0.5]-0.5)/0.5)\n",
    "\n",
    "xlab <- 'temperature'\n",
    "ylab <- 'freq'\n",
    "\n",
    "bins <- seq(0, 50000, 500)\n",
    "\n",
    "# make new colours with transparency\n",
    "or <- alpha('orange', 0.4)\n",
    "bl <- alpha(cbPal[6], 0.4)\n",
    "\n",
    "hist(observed, breaks=bins, xlab=xlab, ylab=ylab, col=bl, main='')\n",
    "hist(sample, breaks=bins, xlab=xlab, ylab=ylab, col=or, add=TRUE)\n",
    "legend('topright', legend=c('sample from theoretical pdf','observed data'), \n",
    "       fill=c(or,bl), cex=0.8, y.intersp=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How can we test more formally whether the observed data appear to be drawn from the theoretical distribution?\n",
    "\n",
    "Here we have a theoretical distribution *which is not directly related to any of the standard statistical distributions*, so parametric methods such as the Shapiro-Wilk test are not applicable.\n",
    "\n",
    "However, we can still perform a *non-parametric* goodness-of-fit test, by comparing the theoretical *cumulative* distributionfunction (cdf) with the empirical one. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "althea_cdf_single <- function(x){\n",
    "    if(x<3000) return(0)\n",
    "    else if(x<4000) return(0.5*(x-3000)/(4000-3000))\n",
    "    else if(x<40000) return(0.5 + 0.5*(x-4000)/(40000-4000))\n",
    "    else return(1)\n",
    "}\n",
    "\n",
    "# We can use the Vectorize() function to make a vector-friendly version of this cdf:\n",
    "althea_cdf <- Vectorize(althea_cdf_single)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The empirical cdf is derived from the observed data using the `ecdf()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "empirical_cdf <- ecdf(observed)\n",
    "\n",
    "x <- seq(0,50000)\n",
    "\n",
    "xlab <- 'temperature'\n",
    "ylab <- 'cumulative probability'\n",
    "    \n",
    "plot(x, althea_cdf(x), type='l', col='orange', xlab=xlab, ylab=ylab)\n",
    "lines(x, empirical_cdf(x), type='l', col=cbPal[6])\n",
    "legend('bottomright', legend=c('theoretical cdf','empirical cdf'), \n",
    "       col=c('orange',cbPal[6]), lty=c(1,1), cex=0.8, y.intersp=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kolmogorov-Smirnov test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Theory\n",
    "\n",
    "The [*Kolmogorov-Smirnov test*](https://en.wikipedia.org/wiki/Kolmogorov–Smirnov_test) (or K-S test) examines the deviation of the empirical cdf, $F_n(x)$, from the theoretical one, $F(x)$, to assess the goodness-of-fit. The test statistic is the quantity\n",
    "\n",
    "$$\n",
    "D_n= \\sup_x |F_n(x)-F(x)|,\n",
    "$$\n",
    "\n",
    "i.e. the greatest vertical distance between the two lines in the plot above.\n",
    "\n",
    "The distribution of $D_n$ under the null hypothesis $H_0: F_n(x) = F(x)$ is called the *Kolmogorov distribution*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Application\n",
    "\n",
    "In our example,\n",
    "\n",
    "$H_0$: $F_n(x) = F(x)$  : The observed temperature distribution of hypergiants is described by Althea's theory.\n",
    "\n",
    "$H_1$: $F_n(x) \\ne F(x)$  : The observed temperature distribution of hypergiants is not described by Althea's theory.\n",
    "\n",
    "$\\alpha = 0.05$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graphically, we can plot $D_n$ for each $x$-value in the plot above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xlab <- 'temperature'\n",
    "ylab <- 'empirical cdf - theoretical cdf'\n",
    "\n",
    "x <- seq(0,50000,1)\n",
    "y_values <- empirical_cdf(x) - althea_cdf(x)\n",
    "\n",
    "plot(x, y_values, type='l', xlab=xlab, ylab=ylab, col='black')\n",
    "\n",
    "D_n <- max(abs(y_values))\n",
    "print(paste('D_n:', D_n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The test itself is easy in R:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result <- ks.test(observed, althea_cdf)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Owing to the rounding in the observed data, there are stars with identical temperatures, i.e. 'ties'. A warning appears because this can affect the calculation of the test statistic. Looking at the cdf plots, it appears unlikely that ties are responsible for much of the deviation between the two lines, so our results should still be reliable.\n",
    "\n",
    "The resulting $p > \\alpha$, so we accept the null hypothesis: the observed data appear to be compatible with Althea's theory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other applications of the K-S test\n",
    "\n",
    "#### Goodness of fit\n",
    "\n",
    "The K-S test can be used as an alternative method for testing normality, or as a goodness-of-fit test for any other theoretical distribution. However, the standard test is only valid if the parameters (e.g. mean and variance) are *not* estimated from the data. \n",
    "\n",
    "If parameters *are* estimated from the data, we will need to use simulation to find an empirical distribution for $D_n$ under $H_0$.\n",
    "\n",
    "\n",
    "#### Two sample test\n",
    "\n",
    "`ks.test()` can also be used to compare two samples (in the absence of a theoretical cdf) to test whether they follow the same distribution. \n",
    "\n",
    "---\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
